!python --version
Dans ce notebook on va essayer de traiter les données des clients d'une banque (ou demandeurs de crédit d'une banque plutôt) avec des outils de Maching Learning, et plus précisément d'apprentissage supervisé.
C'est un apprentissage où les vraies valeurs de la cible nous sont déjà connues. Alors, nos algorithmes sont comme des enfants, ils apprennent d'une partie des données, on valide leur connaissance sur des observations similairas mais pas vues pendant l'apprentissage et ensuite on teste leur performance sur des observations nouvelles. S'ils ont bien appris, mais pas tout appris par coeur sans comprendre (cf. overfitting) on aura des prédictions assez satisfaisantes pour notre cible.
Les deux grandes classes de l'apprentissage supervisé sont:
Classification - cible discrète, en classes
Regression - cible continue
Cet exercice est un exercice de classification. Nous avons des observations qui sont DEFAULT et d'autres qui ne le sont pas, donc il s'agit d'une classification binaire.
Commençons par importer les modules nécessaires pour les tout premières étapes.
# Modules nécessaires
# Essentiels pour le traitement des données
import numpy as np
import pandas as pd
from math import *
from datetime import datetime as dt
# Affichage
import seaborn as sns
import plotly.express as px
from ipywidgets import widgets
from IPython.display import display
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import plotly.figure_factory as ff
pd.options.plotting.backend = "plotly"
# Fixer la graine
np.random.seed(7)
# Ma palette des couleurs
palette = ['#c95c54', '#00bfdc', '#f7ca28', '#39d4a0', '#9c82c7', '#fa9b4d',
'#c95c54', '#00bfdc', '#f7ca28', '#39d4a0', '#9c82c7', '#fa9b4d']
# Fonction pour affichage en pour cent
def to_percent(df):
return df.apply(lambda x: str(round(x*100,2)) + '%')
# Importation des données brutes
df = pd.read_excel("Credit_Cleaned.xlsx")
# Informations sur le dataset
df.info()
# Petit apperçu
df.head()
Nous avons 18 variables explicatives et 1 variable cible.
Variables categorielles:
Les variables qui contiennent des dates:
Discrètes:
Vérifier s'il y a des données non-renseignées. Puisqu'on a la version nettoyée de ces données, on s'attend à ce qu'il n'y ait pas de données nulles, ce qui est bien le cas:
# Potentielles valeurs Null
df.isnull().sum()
Ensuite, sauvegardons les variables categorielles à part:
categorical_features = ['Customer_Type', 'P_Client', 'Source',
'Nb_Of_Products', 'Prod_Sub_Category',
'Educational_Level', 'Marital_Status',
'Type_Of_Residence', 'Number_Of_Dependant',
'Prod_Category', 'Prod_Closed_Date_NA']
descrete_features = ['BirthDate', 'Customer_Open_Date', 'Years_At_Residence',
'Net_Annual_Income', 'Years_At_Business', 'Prod_Decision_Date',
'Prod_Closed_Date']
print(f"Nombre de variables categorielles: {len(categorical_features)}")
print(f"Nombre de variables discrètes: {len(descrete_features)}")
print(f"Au total: {len(categorical_features) + len(descrete_features)} variables explicatives.")
Ces variables doivent être encodées afin qu'on puisse commencer à tester des différents modèles.
Commençons par encoder la cible afin de pouvoir étudier les liens entre elle et les autres variables.
Dans toutes les observations où la valeur de Y est NO_DEFAULT on mettra 0, et là où il y a DEFAULT on mettra 1.
# Encodage de la cible
print(f"Classes avant encodage:\n{df['Y'].value_counts()}\n")
df['Y'].replace({"NO_DEFAULT" : 0, "DEFAULT" : 1}, inplace = True)
print(f"Classes après encodage:\n{df['Y'].value_counts()}")
# Highly skewed data
to_percent(df['Y'].value_counts() / df['Y'].count())
On remarque que la cible est assez déviée (skewed en anglais). C'est-à-dire, on a beaucoup plus d'exemples de NO_DEFAULT que de DEFAULT. Plus précisément, environ 92,7% de NO_DEFAULT et 7.3% de DEFAULT.
Donc, si on prédit 0 pour toute observation d'un ensemble test, on aura une justesse (accuracy) très élevée, mais la précision réelle de l'algorithme sera basse, car on ne va jamais obtenir les prédictions importantes, là où il y a un défault, donc les prédictions de la classe 1.
Ceci va jouer un rôle essentiel dans la définition du modèle qu'on va appliquer à ce jeu de données afin d'obtenir un modèle d'une précision satisfaisante.
D'abord on crée des fonctions qui nous aideront à mieux analyser les variables explicatives (features) et la cible (target).
Commençons par une fonction qui affiche les histogrammes des variables données en argument. On a la possibilité d'indiquer si les variables qu'on veut afficher sont categorielles ou pas grâce au booléan categorical, et aussi un booléan percent qui sert à indiquer si on veut voir les pourcentages que les groupes représentent de toutes les observations. L'argument nb_cols sert à indiquer le nombre de colonnes de l'affichage de tous ces graphiques, à modifier selon les préférances de l'utilisateur.
def plot_features(features, categorical = True, percent = False, nb_cols = 3):
nb_features = len(features)
histnorm_ = ''
if percent: histnorm = 'percent'
# Make subplots
fig = make_subplots(rows=int(np.ceil(nb_features/nb_cols)), cols=nb_cols,
subplot_titles = features)
for i in range(nb_features):
classes = df[features[i]].value_counts()
classes.sort_values(ascending = False, inplace = True)
nb_classes = classes.shape[0]
if nb_classes < 20:
# Counts if categorical feature
fig.add_trace(go.Histogram(x = df[features[i]],
histnorm = histnorm_,
marker_color = palette,
showlegend = False),
row = int(np.floor(i/nb_cols)) + 1, col = i%nb_cols + 1)
else:
# Normal histogram if descrete/continuous feature
fig.add_trace(go.Histogram(x = df[features[i]],
histnorm = histnorm_,
marker_color = palette[4],
showlegend = False),
row = int(np.floor(i/nb_cols)) + 1, col = i%nb_cols + 1)
fig.update_layout(width = 800, height = 800*nb_features/(nb_cols*3))
fig.show()
Affichons donc, toutes les variables, la cible y incluse, pour avoir un premier apperçu de leurs distributions.
plot_features(df.columns,nb_cols = 3)
Ensuite, on créé une fonction qui nous permet d'afficher l'histogramme d'une variable categorielle à gauche et le pourcentage des observations avec Y=0 et Y=1 selon la catégorie dans laquelle elles se trouvent.
def plot_target_vs_categorical_feature(feature):
classes = df[feature].value_counts()
classes.sort_values(ascending = False, inplace = True)
nb_classes = classes.shape[0]
tmp = df.groupby([feature])['Y'].mean()
fig = make_subplots(rows=1, cols=2,
subplot_titles=(f"Client count per {feature}",
"DEFAULT (Y=1) vs. NO_DEFAULT (Y=0) "))
# Count
fig.add_trace(go.Histogram(x = df[feature],
marker_color = palette,
showlegend = False),
row=1, col=1)
# Percentage of target = 1
fig.add_trace(go.Histogram(x = df[feature],
y = df['Y'],
histfunc = 'avg',
# histnorm = 'percent',
name = "Y = 1",
marker_color = palette[0]),
row = 1, col=2)
fig.add_trace(go.Histogram(x = df[feature],
y = 1 - df['Y'],
histfunc = 'avg',
# histnorm = 'percent',
name = 'Y = 0',
marker_color = palette[1]),
row = 1, col=2)
#fig.update_layout(showlegend = False)
fig.update_layout(barmode = 'stack')
fig.update_layout(width = 800, height = 500)
#fig.update_xaxes(categoryorder = 'total descending')
fig.show()
Affichons les 11 variables categorielles et la répartition des valeurs de la cible en fonction des catégories pour chaque variable categorielle.
for feature in categorical_features:
plot_target_vs_categorical_feature(feature)
On dispose de vriament très peu d'observations de DEFAULT. Il semble que le fait d'être un récent locataire (New_Rent) influe sur le fait que Y soit DEFAULT. La même remqrque peut être faite pour la Prod_Category L et pour Prod_Closed_Date_NA = False.
On continue avec les variables discrètes restantes.
# Years_At_Residence
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Residence'],
marker_color = palette[0]))
fig.update_layout(title = "Years_At_Residence")
fig.show()
On doit essayer d'avoir le plus de variables qui ressemblent plutôt à des gaussiennes pour que les modèles (en particulier les modèles linéaires) soient performants. Je décide donc, de modifier cette variable afin d'étaler son histogramme. Je lui applique une transformation avec la fonction racine carrée:
# sqrt(Years_At_Residence)
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Residence'].apply(lambda x: sqrt(x)),
marker_color = palette[6]))
fig.update_layout(title = "Racine carrée de Years_At_Residence")
fig.show()
On voit que les valeurs sont un peu plus centrées et distribuée d'une façon plus normale. Je décide donc de garder cette transformation.
df['Years_At_Residence'] = df['Years_At_Residence'].apply(lambda x: sqrt(x))
# Years_At_Business
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Business'],
marker_color = palette[0]))
fig.update_layout(title = "Years_At_Business")
Cette variable est très déviée également, alors j'applique une transformation racine carrée double. Voici ce que l'on obtient:
# sqrt(sqrt(Years_At_Business))
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Business'].apply(lambda x: sqrt(sqrt(x))),
marker_color = palette[0]))
fig.update_layout(title = "Years_At_Business")
On la garde également.
df['Years_At_Business'] = df['Years_At_Business'].apply(lambda x: sqrt(sqrt(x)))
# Net Annual Income
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Net_Annual_Income'],
marker_color = palette[0]))
Cette variable est extrêmement déviée, très probablement parce qu'il y a très peu de gens qui ont des gros salaires, comparé au nombre de gens qui sont pauvres, ou payés correctement. ~Vive le capitalisme~.
Alors je la modifie en lui appliquant une transformation log.
# log(Net Annual Income)
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Net_Annual_Income'].apply(lambda x: log(x,10)),
marker_color = palette[0]))
On décide donc, de prendre plutôt le logarithme en base 10 du salaire net annuel pour nos modèles éventuels.
df['Net_Annual_Income'] = df['Net_Annual_Income'].apply(lambda x: log(x,10))
plot_features(['BirthDate',
'Customer_Open_Date',
'Prod_Decision_Date',
'Prod_Closed_Date'], nb_cols=1)
Créons une nouvelle variable qui va s'appeler duration qui contiendra la période entre Prod_Closed_Date et Customer_Open_Date. Ainsi, on espère se libérer des données très déviées tout en gardant l'information importante qu'elles portent.
duration = df['Prod_Closed_Date'] - df['Customer_Open_Date']
duration = duration.dt.days
duration.name = 'days'
Regardons son histogramme:
ig = go.Figure()
fig.add_trace(go.Histogram(x = duration,
marker_color = palette[0],
showlegend = False))
On remarque qu'elle est assez déviée, alors on applique une transformation log.
fig = go.Figure()
fig.add_trace(go.Histogram(x = duration.apply(lambda x: log(x,10)),
marker_color = palette[0],
showlegend = False))
On décide d'enlever les variables Customer_Open_Date et Prod_Decision_Date et remplacer l'information qu'elles apportent par leur différence, une durée, plus précisément on prend le logarithme en base 10 de ces durées pour que ça ressebmle le plus à une loi gaussienne.
duration = duration.apply(lambda x: log(x,10))
df.drop(columns = ['Customer_Open_Date', 'Prod_Closed_Date'], inplace = True)
df['prod_duration'] = duration
Créons une fonction qui nous permet de convertir les données de format datetime en un format ordinal.
def convert_datetime_featues_to_ordinals(df, features):
for feature in features:
tmp = df[feature].apply(lambda x: x.toordinal())
tmp = (tmp - tmp.mean()) / tmp.std()
df[feature] = tmp.values
convert_datetime_featues_to_ordinals(df, ['BirthDate', 'Prod_Decision_Date'])
df.columns
Voici un petit aperçu de notre data set à présent:
df.head()
corr = df.corr()
mask = np.triu(np.ones_like(corr, dtype=bool))
px.imshow(corr, color_continuous_scale = 'viridis')
Pour cela, on va encoder les variables Customer_Type et Prod_Closed_Date à valeurs dans {0,1} et on va utiliser la méthode OneHotEncoding là où c'est nécessaire.
Classes pour lesquelles il faudrait encoder les données:
Variables à 2 classes:
Variables à 3 classes:
Variables à 4 ou plus classes:
# Nombre de classes pour les variables categorielles
for column in df.columns:
classes = df[column].value_counts()
if classes.count() < 20:
print(f"Variable: {column}\nNombre de classes: {classes.count()}\nClasses:\n{classes}\n\n")
# Encodage de Customer_Type
print(f"Classes avant encodage:\n{df['Customer_Type'].value_counts()}\n")
df['Customer_Type'].replace({"Non Existing Client" : 0, "Existing Client" : 1}, inplace = True)
print(f"Classes après encodage:\n{df['Customer_Type'].value_counts()}")
# Encodage de P_Client
print(f"Classes avant encodage:\n{df['P_Client'].value_counts()}\n")
df['P_Client'].replace({'NP_Client' : 0, 'P_Client' : 1}, inplace = True)
print(f"Classes après encodage:\n{df['P_Client'].value_counts()}")
# Encodage de Prod_Closed_Date_NA
print(f"Classes avant encodage:\n{df['Prod_Closed_Date_NA'].value_counts()}\n")
df["Prod_Closed_Date_NA"].replace({True : 1, False : 0}, inplace = True)
print(f"Classes après encodage:\n{df['Prod_Closed_Date_NA'].value_counts()}")
# Encodage de Source
print(f"Classes avant encodage:\n{df['Source'].value_counts()}\n")
df["Source"].replace({'Sales' : 1, 'Branch' : 0}, inplace = True)
print(f"Classes après encodage:\n{df['Source'].value_counts()}")
Il nous reste à encoder les variables categorielles restantes, à savoir:
def OneHotEncoding_(df, features):
for feature in features:
tmp = pd.get_dummies(df[feature])
for column in tmp.columns:
df[column] = tmp[column].values
df.drop(columns = features, inplace = True)
features = ['Educational_Level', 'Marital_Status', 'Prod_Sub_Category',
'Type_Of_Residence', 'Prod_Category']
OneHotEncoding_(df, features)
Une fois nos variables bien choisies, modifiées et encodées, on peut procéder à l'étape suivante et choisir le modèle de prédiction qui conviendrait le mieux.
df.columns
df.dtypes
Tout d'abord, séparons notre dataset en features (X) et target (y)
Remarque: Je prends des copies car la dataframe est un objet mutable en python, donc il se comporte comme un passage par référence dans un langage compilé. Puisque je ne veux pas modifier la dataframe initiale, je prends des copies.
X = df.drop(columns = ['Y']).copy()
y = df['Y'].copy()
Ensuite, on va partager notre dataset en trois parties, une sur laquellue on va entraîner (X_train), une qui va nous servir pour la cross-validation (X_test) et une troisième qui va nous servir pour tester les derniers résultats (X_new_data). On va également standardiser les features pour nous assurer que les algorithmes de descente de gradient optimisés (qui se cachent derrière presque tous les modèles de machine learning) puissent converger.
Pour cela, on va se servir d'une utilité du module sklearn:
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.45)
X_test, X_new_data, y_test, y_new_data = train_test_split(X_test, y_test, test_size=0.25)
sc = StandardScaler()
sc.fit(X_train)
X_train = sc.transform(X_train)
X_test = sc.transform(X_test)
X_new_data = sc.transform(X_new_data)
from sklearn import metrics
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
model = LogisticRegression(random_state=0)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
Lorsque la cible est si déviée, il est très important qu'on observe la matrice de confusion, qui est de la forme:
vraies négatifs | faux positifs
faux négatifs | vraies positifs
Les vraies négatifs sont les personnes qui ne présentent pas de défault sur leur crédit et dont on a également prédit qu'il n'y avait pas de défaut.
Les vraies positifs sont les personnes qui présentent de défault sur leur crédit et dont on a également prédit qu'il y avait de défault.
Les faux négatifs sont les personnes pour lesquelles on avait prédit qu'elles ne présenteraient pas de défaut, mais pourtant elles ont le défault.
Les faux_positifs sont les personnes qui ne présentent pas de défault, mais pour lesquelles on avait prédit qu'elles présentaient de défault, à tort.
Regardons ces valeurs pour ce modèle.
true_positives = (((y_test == 1) & (y_pred == 1)) == True).sum()
false_positives = (((y_test == 0) & (y_pred == 1)) == True).sum()
true_negatives = (((y_test == 0) & (y_pred == 0)) == True).sum()
false_negatives = (((y_test == 1) & (y_pred == 0)) == True).sum()
print(f"Vrais positifs: {true_positives}")
print(f"Faux positifs: {false_positives}")
print(f"Vrais négatifs: {true_negatives}")
print(f"Faux négatifs: {false_negatives}")
Il existe une fonction dans le module sklearn, dans la classe metrics, qui nous donne la matrice de confusion directement.
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(y_test, y_pred)
print(cm)
Il est pratique de pouvoir bien visualiser cette matrice de confusion
def plot_confusion_matrix(cm):
z_text = [[f'True Negatives = {cm[0,0]}',
f'False Positives = {cm[0,1]}'],
[f'False Negatives = {cm[1,0]}',
f'True Positives = {cm[1,1]}']]
fig = ff.create_annotated_heatmap(cm,
x=['0','1'], y=['0','1'],
annotation_text=z_text, colorscale='plasma')
# add custom xaxis title
fig.add_annotation(dict(font=dict(color="black",size=15),
x=0.5,
y=-0.15,
showarrow=False,
text="Predicted value",
xref="paper",
yref="paper"))
# add custom yaxis title
fig.add_annotation(dict(font=dict(color="black",size=15),
x=-0.12,
y=0.5,
showarrow=False,
text="Real value",
textangle=-90,
xref="paper",
yref="paper"))
# adjust margins to make room for yaxis title
fig.update_layout(margin=dict(t=50, l=130))
# add colorbar
fig['data'][0]['showscale'] = True
fig.show()
plot_confusion_matrix(cm)
Pour bien étudier la qualité de ce modèle il est important de calculer la precision et le recall. Ces quantités sont définies ainsi:
PRECISION : De toutes les observations pour lesquelles on a prédit un défault (y_pred = 1), quelle partie a vraiment un défaut (y_test = 1)?
RECALL : Pour quelle partie des observations qui ont un défaut a-t-on réussi à prédire le défaut?
Il est assez intuitif à voir que si on prédit 0 pour toute observation du test on aura une précision à 100%, et une justesse de 92.3% (le pourcentage des observations qui ne présentent pas de défault) mais ce modèle n'a aucun intérêt général.
De même, si on prédit 1 pour toute observation, on aura un recall à 100%. Mais, idem, ce modèle est évidemment inutil.
Donc, notre modèle sera bon s'il a une grande précision et un grand recall.
C'est pour cela qu'on utilise le F-score, également appelé $F_1$. En nottant P la précision et R le recall, le $F_1$ est donné par:
Créons une fonction donc, qui nous calculera ces valeurs pour nos données et nos modèles.
def precision_recall_fscore_accuracy(y_test, y_pred):
true_positives = (((y_test == 1) & (y_pred == 1)) == True).sum()
false_positives = (((y_test == 0) & (y_pred == 1)) == True).sum()
true_negatives = (((y_test == 0) & (y_pred == 0)) == True).sum()
false_negatives = (((y_test == 1) & (y_pred == 0)) == True).sum()
print(f"Vrais positifs: {true_positives}")
print(f"Faux positifs: {false_positives}")
print(f"Vrais négatifs: {true_negatives}")
print(f"Faux négatifs: {false_negatives}")
precision = true_positives / (true_positives + false_positives)
recall = true_positives / (true_positives + false_negatives)
print(f"Precision: {precision:.3f}\nRecall: {recall:.3f}")
F1_score = 2 * (precision * recall) / (precision + recall)
print(f"F1 Score: {F1_score:.3f}")
accuracy = (true_positives + true_negatives) / y_test.shape[0]
print(f"Accuracy: {accuracy:.3f}")
return precision, recall, F1_score, accuracy
precision_recall_fscore_accuracy(y_test, y_pred);
Pour augmenter la probabilité de détecter les défauts on va changer le seuil à partir duquel on décide si on prédit 0 ou 1 dans la fonction sigmoid dans la régression logistique. Pour cela, on créé un vecteur de différents seuils et on compare le $F_1$ score pour chacun. A la fin, on choisit celui qui maximise le $F_1$.
thresholds_ = np.linspace(0,1,30)[1:-1]
F1_score_list = []
for threshold in thresholds_:
y_pred_ = (model.predict_proba(X_test)[:,1] >= threshold).astype(bool)
precision, recall, F1_score, accuracy = precision_recall_fscore_accuracy(y_test, y_pred_);
F1_score_list.append(F1_score)
tmp = pd.DataFrame(F1_score_list, index = thresholds_, columns = ['F1_score'])
tmp.dropna(inplace = True)
best_threshold_ = tmp.iloc[np.argmax(tmp)]
best_threshold_
fig = go.Figure()
fig.add_trace(go.Scatter(x = thresholds_, y = F1_score_list,
marker_color = palette[0]))
fig.update_layout(title = "F1 Score en fonction de différents seuils")
fig.update_xaxes(title = "Threshold")
fig.update_yaxes(title = "F1-Score")
On obtient que le meilleur seuil est environ 0.345. Appliquons cela à nos prédictions
print("Threshold = 0.5")
precision_recall_fscore_accuracy(y_test, y_pred);
print("\nThreshold = 0.38")
y_modified_pred = (model.predict_proba(X_test)[:,1] >= 0.345).astype(int)
precision_recall_fscore_accuracy(y_test, y_modified_pred);
Ceci n'augmente pas le $F_1$ de façon importante, mais si on veut plutôt avoir meilleure probabilité de reconnaître le plus de vrais défaults possibles (au prix des faux positifs bien évidemment) alors cela serait la stratégie à adopter.
Ensuite traçons la courbe d'apprentissage, calculée sur 10 validations croisées, en utilisant tous les processeurs (parallélisation), qui montre la valeur moyenne du F score sur le train, et sur le test, en fonction de la taille du train.
# Learning Curve
from sklearn.model_selection import learning_curve
train_sizes, train_scores, test_scores = learning_curve(LogisticRegression(),
X, y,
cv = 10,
scoring = 'f1',
n_jobs=-1,
train_sizes=np.linspace(0.01, 1.0, 50))
train_means = np.mean(train_scores, axis=1)
test_means = np.mean(test_scores, axis=1)
fig = go.Figure()
fig.add_trace(go.Scatter(x = train_sizes, y = train_means,
marker_color = palette[1],
name = 'J_train'))
fig.add_trace(go.Scatter(x = train_sizes, y = test_means,
marker_color = palette[0],
name = 'J_cv'))
fig.update_xaxes(title = 'Train size')
fig.update_yaxes(title = 'F_score')
fig.update_layout(title = "Learning curve")
Il semble qu'on a un grand biais, mais ceci est attendu vu que notre cible est très biaisée.
train_sizes, train_scores, test_scores = learning_curve(LogisticRegression(),
X, y,
cv = 10,
scoring = 'roc_auc',
n_jobs=-1,
train_sizes=np.linspace(0.01, 1.0, 50))
train_means = np.mean(train_scores, axis=1)
test_means = np.mean(test_scores, axis=1)
fig = go.Figure()
fig.add_trace(go.Scatter(x = train_sizes, y = train_means,
marker_color = palette[1],
name = 'Train'))
fig.add_trace(go.Scatter(x = train_sizes, y = test_means,
marker_color = palette[0],
name = 'CV'))
fig.update_xaxes(title = 'Train size')
fig.update_yaxes(title = 'ROC_AUC')
fig.update_layout(title = "Learning curve")
Pour remédier au probllème de grand biais, on pourrait faire les choses suivantes:
On essayera de voir comment on pourrait augmenter C de façon intelligente. Voici une routine qui teste plusieurs valeurs différentes de C et qui calcule le $F_1$ pour chaque modèle appliqué avec chaque valeur de C. Ensuite, on n'aura qu'à choisir le paramètre C qui maximise le $F_1$
# Choix du meilleur paramètre de régularisation C
from sklearn.metrics import f1_score
C_values = np.linspace(0.05,12,30)
f1_scores_train = []
f1_scores_cv = []
for C in C_values:
model = LogisticRegression(C = C, random_state=0).fit(X_train, y_train)
y_pred = model.predict(X_test)
y_pred_train = model.predict(X_train)
f1_scores_cv.append(f1_score(y_test,y_pred, zero_division=1))
f1_scores_train.append(f1_score(y_train,y_pred_train, zero_division=1))
fig = go.Figure()
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_train,
marker_color = palette[1],
name = 'Train'))
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_cv,
marker_color = palette[0],
name = 'CV'))
fig.update_xaxes(title = 'Regularization parameter C')
fig.update_yaxes(title = 'Scores')
fig.update_layout(title = "Learning curve")
Malheureusement, on remarque que la valeur de C n'a pas de grande influence sur le $F_1$, la valeur par défaut (qui est égale à 1) est donc suffisante et on ne changera pas ce paramètre.
Affichons la courbe ROC qui représente le Taux de Vrais Positifs (TVP) (ce qui n'est rien d'autre qu'un synonime pour recall) en fonction du Taux de Faux Positifs (TFP), calculée à partir de la fonction décision et les vraies labels de la cible.
Soit $TVP$ le Taux de Vrais Positifs, $TFP$ le Taux de Faux Positifs, et donc $V$ ou $F$ signifierait Vrais ou Faux, et $P$ et $N$ positifs et négatifs, respectivement. Nous avons:
$$ TVP = \frac{VP}{VP + FN}$$$$ TFP = \frac{FP}{FP + VN}$$Affichons également la courbe Précision/Recall qui représente les valeurs de la précision en fonction des valeurs du recall (selon les seuils: thresholds). L'intégrale de cette courbe est en fait une très bonne mesure de performance pour les modèles de classification: si toutes les prédictions du modèles sont fausses cette valeur sera égale à 0, si toutes les prédictions sont bonnes elle sera égale à 1. Néanmoins, ici il faut se méfier car la cible est très fortement biaisée. L'AUC est invariante aux changements de seuil de choix (threshold). Donc, il faut se méfier dans les cas où il est préférable de maximiser les vrais positifs, quitte à avoir plus de faux positifs (comme cet exercice où il est préférable de trouver tous les défaults potentiels) ou par exemple lorsqu'on veut minimiser le nombre de faux positifs, quitte à avoir plus de faux négatifs (dans le cas de données médicales et de maladies terminales par exemple).
# ROC and Precision/Recall Curve
from sklearn.metrics import roc_curve, auc, precision_recall_curve, average_precision_score
y_score = model.fit(X_train,y_train).decision_function(X_test)
false_positive_rate, true_positive_rate, _ = roc_curve(y_test, y_score)
roc_auc = auc(false_positive_rate, true_positive_rate)
ap = average_precision_score(y_test, y_score)
precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fig = make_subplots(rows = 1, cols = 2,
subplot_titles = [f"ROC curve (AUC = {roc_auc:.2f})",
f"Precision-Recall curve (AP = {ap:.2f})"])
fig.add_trace(go.Scatter(x = false_positive_rate, y = true_positive_rate,
marker_color = palette[0]),
row = 1, col = 1)
fig.add_trace(go.Scatter(x = recall, y = precision,
marker_color = palette[1]),
row = 1, col = 2)
fig.update_layout(showlegend = False)
Les résultats sont satisfaisants, je choisirais quand même de garder la modification du seuil pour maximiser le nombre de défaults trouvés. Le modèle de régression logistique que je garde serait donc:
logr = LogisticRegression(random_state = 0, C = 1).fit(X_train, y_train)
On va appliquer le modèle SVM pour la classification, donc plus connu sous le nom de SVC
model = SVC(kernel = 'poly', random_state = 3)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
precision_recall_fscore_accuracy(y_test, y_pred);
y_score = model.fit(X_train,y_train).decision_function(X_test)
false_positive_rate, true_positive_rate, _ = roc_curve(y_test.values, y_score)
roc_auc = auc(false_positive_rate, true_positive_rate)
ap = average_precision_score(y_test, y_score)
precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fig = make_subplots(rows = 1, cols = 2,
subplot_titles = [f"ROC curve (AUC = {roc_auc:.2f})",
f"Precision-Recall curve (AP = {ap:.2f})"])
fig.add_trace(go.Scatter(x = false_positive_rate, y = true_positive_rate,
marker_color = palette[0]),
row = 1, col = 1)
fig.add_trace(go.Scatter(x = recall, y = precision,
marker_color = palette[1]),
row = 1, col = 2)
fig.update_layout(showlegend = False)
Ce modèle se comporte moins bien que la régression logistique, j'avais testé tous les kerneles possibles et celui là donnait les meilleurs résultats, cependant pas assez suffisants. Essayons de modifier le paramètre de régularization.
C_values = np.linspace(0.05,30,50)
from sklearn.metrics import roc_auc_score, f1_score
f1_scores_train = []
f1_scores_cv = []
for C in C_values:
model = SVC(C = C, kernel = 'poly', random_state = 3).fit(X_train, y_train)
y_pred = model.predict(X_test)
y_pred_train = model.predict(X_train)
f1_scores_cv.append(f1_score(y_test,y_pred, zero_division=1))
f1_scores_train.append(f1_score(y_train,y_pred_train, zero_division=1))
fig = go.Figure()
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_train,
marker_color = palette[1],
name = 'Train'))
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_cv,
marker_color = palette[0],
name = 'CV'))
fig.update_xaxes(title = 'Regularization parameter C')
fig.update_yaxes(title = 'F Score')
fig.update_layout(title = "F Score in function of Regularization parametar C")
Ici on remarque une amélioration. Plus on régularise, mieux c'est, avec une stabilisation autour de 14. Donc, on choisit la valeur 14 pour C.
model = SVC(kernel = 'poly', C = 14)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
precision_recall_fscore_accuracy(y_test, y_pred);
Observons la matrice de confusion:
cm = confusion_matrix(y_test, y_pred)
plot_confusion_matrix(cm)
Et finalement, affichons la courbe ROC encore une fois:
y_score = model.decision_function(X_test)
false_positive_rate, true_positive_rate, _ = roc_curve(y_test.values, y_score)
roc_auc = auc(false_positive_rate, true_positive_rate)
ap = average_precision_score(y_test, y_score)
precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fig = make_subplots(rows = 1, cols = 2,
subplot_titles = [f"ROC curve (AUC = {roc_auc:.2f})",
f"Precision-Recall curve (AP = {ap:.2f})"])
fig.add_trace(go.Scatter(x = false_positive_rate, y = true_positive_rate,
marker_color = palette[0]),
row = 1, col = 1)
fig.add_trace(go.Scatter(x = recall, y = precision,
marker_color = palette[1]),
row = 1, col = 2)
fig.update_layout(showlegend = False)
On se rend compte que selon ces métriques là, la performance a baissé. Mais, le recall a augmenté, et avec cela le $F_1$ aussi. Donc, on reste dans le cadre d'un modèle satisfaisant à nos attentes. On le garde dans:
svc = SVC(kernel = "poly", C = 14, random_state = 3).fit(X_train, y_train)
Le troisième, et dernier modèle qu'on va essayer c'est le modèle de Forêt Aléatoire.
Les hypermarapètres ici, je les ai cherché intuitivement car j'ai de l'expérience avec des données similaires, autant en quantité qu'en qualité, et le plus important - en biais. Voici ce que l'on obtient:
from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier(n_estimators = 150,
max_depth = 13,
bootstrap = False,
min_samples_split = 8,
random_state = 7,
n_jobs=-1)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
precision_recall_fscore_accuracy(y_test, y_pred);
cm = confusion_matrix(y_test, y_pred)
plot_confusion_matrix(cm)
Ce modèle a de très bonnes pérformances sur ces données. Comparable à la régression logistique qui était très pérformante également.
from sklearn.metrics import classification_report
print(classification_report(y_test, y_pred))
from sklearn.metrics import roc_auc_score
roc_auc_score(y_test, y_pred)
Le AUC est plus petit que celui des autres modèles, mais comme mentionné précédemment, cela n'est pas une métrique complètement bien adapté à notre cas.
On sauvegarde ce modèle et on procède à la dernière étape.
rf = RandomForestClassifier(n_estimators = 150,
max_depth = 13,
bootstrap = False,
min_samples_split = 8,
random_state = 7,
n_jobs=-1).fit(X_train, y_train)
Voici une comparaison des modèles sur les données test qu'on a séparées avant.
# Comparaison des modèles
print("---------------- Logistic Regression -----------------")
logr_pred = logr.predict(X_new_data)
print(classification_report(y_new_data, logr_pred))
print(f"AUC = {roc_auc_score(y_new_data, logr_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, logr_pred))
print("----------------------- SVC ---------------------------")
svc_pred = svc.predict(X_new_data)
print(classification_report(y_new_data, svc_pred))
print(f"AUC = {roc_auc_score(y_new_data, svc_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, svc_pred))
print("------------------- Random Forest ---------------------")
rf_pred = rf.predict(X_new_data)
print(classification_report(y_new_data, rf_pred))
print(f"AUC = {roc_auc_score(y_new_data, rf_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, rf_pred))
On peut conclure que d'après plusieurs critères de performance, le modèle de Régression Logistique est le plus performant pour ces données.
Afin de privilégier la détéction des défaults, appliquons la baisse du seuil comme précédamment vu:
logr_modified_pred = (logr.predict_proba(X_new_data)[:,1] >= 0.345).astype(int)
precision_recall_fscore_accuracy(y_new_data, logr_modified_pred)
print(f"AUC = {roc_auc_score(y_new_data, logr_modified_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, logr_modified_pred))
Il existe plusieurs modèles différentes pour l'apprentissage supervisé, et peut importe celui qu'on choisi, l'important c'est de pouvoir bien intérpréter les données. Même si on n'est pas familier avec toutes les variables explicatives en détails, il faudrait toujours garder une intuition réveillée en permanance, avant de ne pas se perdre dans les million d'options et paramètres qui existent et de garder un esprit critique tout au long de l'anaylse, la modélisation et la validation.
Pour cet exercice, il existe plusieurs façon d'améliorer le modèle de régression logistique ou les autres modèles. On pourrait choisir une méthode différente pour la régularisation, on pourrait construire de nouvelles variables explicatives polynômiales à partir de celles dont on dispose déjà, on pourrait lancer un RandomGridSearch en parallèle, avec une K-fold validation croisée... Cela peut être le travail de plusieurs mois si on court vraiment derrièrre ne serait ce que quelques petits pourcents de précision (accuracy) de plus.